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In this paper we derive one- and two-sample multivariate em- 
pirical Bayes statistics (the MB-statistics) to rank genes in order of 
interest from longitudinal replicated developmental microarray time 
course experiments. We first use conjugate priors to develop our 
one-sample multivariate empirical Bayes framework for the null hy- 
pothesis that the expected temporal profile stays at 0. This leads to 
our one-sample MB-statistic and a one-sample T 2 -statistic, a vari- 
ant of the one-sample Hotelling T 2 -statistic. Both the MB-statistic 
and T 2 -statistic can be used to rank genes in the order of evidence 
of nonzero mean, incorporating the correlation structure across time 
points, moderation and replication. We also derive the corresponding 
MB-statistics and T 2 -statistics for the one-sample problem where the 
null hypothesis states that the expected temporal profile is constant, 
and for the two-sample problem where the null hypothesis is that two 
expected temporal profiles are the same. 



1. Introduction. Microarray time course experiments differ from other 
microarray experiments in that gene expression values at different time 
points can be correlated. This may happen when the design is longitudi- 
nal, that is, where the mRNA samples at successive time points are taken 
from the same units. Such longitudinal experiments make it possible to mon- 
itor and study the temporal changes within units of biological processes of 
interest for thousands of genes simultaneously. Two major categories of time 
course experiments are those involving periodic and developmental phenom- 
ena. Periodic time courses typically concern natural biological processes such 
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as cell cycles or circadian rhythms, where the temporal profiles follow reg- 
ular patterns [7, 8, 25, 26]. On the other hand, developmental time course 
experiments measure gene expression levels at a series of times in a de- 
velopmental process, or after applying a treatment such as a drug to the 
organism, tissue or cells [9, 32, 34]. In this case, we typically have few prior 
expectations concerning the temporal patterns of gene expression. The gene 
ranking methods we develop in this paper are mainly for longitudinal repli- 
cated developmental time course data. 

A typical microarray time course dataset consists of expression measure- 
ments of G genes across k time points, under one or more biological condi- 
tions (e.g., wildtype versus mutant). The number of genes G (1,000-40,000) 
is very much larger than the number of time points k, which can be 5-10 
for shorter and 11-20 for longer time courses. Many such experiments are 
unreplicated due to cost or other limitations, and when replicates are done, 
the number n is typically quite small, say, 2-5. We refer the reader to [29] 
for a fuller review of microarray time course experiments. 

One of the statistical challenges here is to identify genes of interest. In 
what we call the one-sample problem, these are genes whose patterns of ex- 
pression change over time, perhaps in some specific way. In the two-sample 
problem we seek genes whose temporal patterns differ across two biological 
conditions. Such genes are of interest to biologists because they are often in- 
volved in the biological processes motivating the experiment. The challenge 
arises from the fact that there are very few time points, and very few repli- 
cates per gene. The series are usually so short that we cannot consider using 
standard time series methods as described in [10], such as Fourier analysis, 
ARMA models or wavelets. The methods proposed in this paper are for the 
one- and two-sample problems with longitudinal replicated microarray time 
course experiments of the developmental kind. 

The gene ranking problem for such microarray experiments is relatively 
new. Few methods have been proposed to deal specifically with these prob- 
lems. The most widely used method for identifying temporally changing 
genes in replicated microarray experiments is to carry out multiple pairwise 
comparisons across times, using statistics developed for comparing two inde- 
pendent samples, [2, 6, 13, 19, 20, 23, 24, 33]. These methods are not entirely 
appropriate as they do not incorporate the fact that longitudinal microar- 
ray time course samples are correlated. A simple and intuitive approach 
to our problem is to use classical or mixed ANOVA models; see Chapter 
6 of [11] for a discussion of the latter for analyzing longitudinal data and 
[22] for a modified approach based on the former for use in the microarray 
context. However, a number of questions are not adequately addressed by 
the classical ANOVA methods, or the variants of [22]. As with the pairwise 
comparisons, the -F-statistic assumes that gene expression measurements 
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at different times are independent. The classical ANOVA models also as- 
sume normality of the gene expression measurements, which may not be a 
great concern when these are on the log scale. More importantly, standard 
.F-statistics in this context will generally lead to more false positives and 
false negatives than is desirable, due to poorly estimated error variances in 
the denominator. This issue can be dealt with using the idea of moderation; 
see, for example, [2, 14, 20, 24, 33]. Moderation in our longitudinal context 
means the smoothing of gene-specific sample variance-covariance matrices 
toward a common matrix. When we do this, fewer genes which are not dif- 
ferentially expressed over time but have very small replicate variances are 
falsely identified as being differentially expressed, and fewer genes which are 
differentially expressed over time but have large replicate variances (e.g., 
due to outliers) are missed by the F-statistic. 

A moderated gene-specific score based on the Wald statistic for the longi- 
tudinal one-sample problem was proposed in [15]. However, their method is 
only applicable to the situation when the number of replicates is greater than 
the number of time points. In [3] the expression profiles for each gene and 
each of two biological conditions were represented by continuous curves fitted 
using B-splines. A global difference between the two continuous curves and 
an ad hoc likelihood-based p-value was calculated for each gene. B-splines 
were also used in [17] to identify genes with different temporal profiles in the 
two-sample case. Recently, B-splines were again adopted by [27] to model 
the population mean, constructing the ^-statistics for both longitudinal and 
cross-sectional data with one or more biological conditions. A major feature 
of this paper was a careful treatment of the multiple testing issue in this 
context. A novel HMM approach which incorporates the dependency in gene 
expression measurements across times was proposed in [36] for data with two 
or more biological conditions. This is one example of using HMM to identify 
differentially expressed genes across at least two biological conditions in this 
context. 

The multivariate empirical Bayes model proposed in this paper was mo- 
tivated by the analogous univariate model proposed in [20] for identifying 
differentially expressed genes in two-color comparative microarray experi- 
ments, and the more recent extensions by Smyth [24]. The 5-statistic in 
[20] and [24], and the univariate moderates £-statistic t g in [24], consider 
just one parameter or contrast at a time in the null hypotheses. They are 
not for null hypotheses with two or more parameters or contrasts of inter- 
est simultaneously. However, a partly-moderated ^-statistic was introduced 
in [24], which moderates the error variance in the denominator of the or- 
dinary F. This partly-moderated ^-statistic is useful for the simultaneous 
comparison of multiple uncorrelated contrasts, but as mentioned above, it 
is not appropriate for longitudinal experiments. Both the MB-statistics and 
the T 2 -statistics derived in this paper provide a degree of moderation, while 
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retaining the temporal correlation structure. They can be used with both 
single- and two-channel microarray experiments. 

This paper is organized as follows. After briefly introducing our notation, 
we formally state the null and alternative hypotheses for the gene ranking 
problem in Section 2. In Section 3 we present moderated versions of the stan- 
dard likelihood-ratio and Hotelling T 2 -statistics. We formally build up our 
multivariate empirical Bayes model and derive the MB- and T 2 -statistics 
in Section 4. A brief description of a case study is presented in Section 5. 
Section 6 reports results from a simulation study in which we compare the 
one-sample MB-statistic (the T 2 -statistic) of Section 4.3 with other statis- 
tics. We discuss our models and give directions for future work in Section 7. 

Now we introduce some notation for one-sample problems. For two-sample 
problems, the notation is similar and easily conceived. For each gene g, 
g = 1, . . . ,G, we have n g independent time series, and we model these as 
i.i.d. kxl random vectors ~K g i, . . . ,~K grig with gene-specific means fi g and 
covariance matrices S 9 . Since relative or absolute gene expression measure- 
ments are approximately normal on the log scale, we make the multivariate 
normality assumption on data X 5 i , . . . , X 5 „, 9 . Our results are to be judged on 
their practical usefulness, not on the precise fit of our data to a multivariate 
normal distribution. As will be seen shortly, our final formulae involve the 
multivariate t distribution. Thus, a measure of robustness is built in, and 
our approach will probably be about as effective for elliptically distributed 
random vectors. We use the natural conjugate priors for n g and that is, 
an inverse Wishart prior for H g and a dependent multivariate normal prior 
for n g . To simplify the notation, the subscript g will be dropped for the rest 
of this paper. The statistical models presented in the rest of this paper are 
for an arbitrary single gene g. 

The details in this paper differ in two ways from the standard conjugate 
priors. First, we also have an indicator / such that 1=1 when the alternative 
K is true and 1 = when the null H is true, with the priors for n differing in 
these two cases. Second, when the null hypothesis states that a gene's mean 
expression level is constant, in order to get a simple closed form expression 
for the posterior odds, we assume that the gene-specific covariance matrix H 
commutes with the k x k projection matrix P = A; _1 lfcl^,, that is, PS = SP. 
In this case the k x k inverse Wishart prior for 5] is replaced by a (k — 1) x 
(k — 1) inverse Wishart prior for a part of S and an inverse gamma prior 
for the remainder. These two-part priors are independent; see Section 4.3 
for details. 

2. Hypothesis testing. Our gene ranking problem will be formally stated 
as a hypothesis testing problem. In this paper we only seek a statistic for 
ranking genes in the order of evidence against the null hypothesis; we do 
not hope to obtain raw or adjusted p- values as in [27]. 
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Following the notation in [4], the null hypothesis is denoted by H, while 
the alternative hypothesis is denoted by K. The null hypothesis correspond- 
ing to a gene's mean expression level being is H : \x = 0, 5] > 0; the alter- 
native is IT : /x 7^ 0, £ > 0. An easy extension is H : fi = fj, , X > versus 
K : [i ^ fj, , S > 0, where ytt is known. Later we will consider the null hy- 
pothesis that a gene's expression is constant, against the alternative that it 
is not: H :^l = S > 0, where fio is a scalar representing the expected 
value of the gene's expression level at any time point under H, and 1 is 
the k x 1 constant vector of Is; K : [j,^ (iqI, S > 0. Finally, we will consider 
the null hypothesis that a gene's mean expression levels are the same under 
two different biological conditions, versus the alternative that they are not: 
H:n z = fi Y , V Z = -£ Y = £>0; K:fi z ^n Y , £ z = £y = £ > 0. 

3. The moderated Xi?-statistic. 

3.1. One-sample or paired two-sample problem. A likelihood-ratio statis- 
tic can be used directly to test the null hypothesis H against the alternative 
hypothesis K when n> k. According to standard multivariate results (e.g., 
[21]), under the alternative hypothesis that there are no constraints on /j, 
and XI, their maximum likelihood estimates are 

A a- = X, 

Tik = S, 

n 

where S = (n — 1) _1 Y^}=i 0^4 ~ X) (Xj — X)' is the sample variance-covariance 
matrix. Also as in [21], the maximum likelihood estimate for the uncon- 
strained X! under the null that fi = fi is 

n 

where d = (i K — fi H , the difference between the maximum likelihood esti- 
mates for fj, (if unknown) under H and K. If fi is known, then d = fj, K — fj, . 
The likelihood ratio statistic for testing any such H against the above K is 

LR = 2(J£ ax - fg 8 *) 

(3-1) 

= nlog(l + -^d'S _1 d 

If the null hypothesis states that fi = 0, then d reduces to X, hence, the 
likelihood-ratio statistic is equation (3.1) with d replaced by X. Similarly, 
if fj, is known, then d = X — fi . On the other hand, if the null hypothesis 
states that /j, = [J,ol, then the maximum likelihood estimate for /i is 

/l'S" 1 ^ 
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The statistic nd'S~ 1 d is the one-sample Hotelling T 2 -statistic, and by Sec- 
tion 5.3.1b in [21], it is distributed as a Hotelling T 2 (k — l,n — 1) under H, 
that is, ((n — k + l)nd'S _1 d)/((n — l)(fc — 1)) has an F-distribution with 
degrees of freedom [k — l,n — k + 1). 

In the microarray time course context, the number of replicates n is typ- 
ically smaller than the number of time points k, and so S has less than full 
rank. Furthermore, as discussed in [29], we wish to moderate the sample 
variance-covariance matrix. Our moderated S will take the form 

~_ wl + (ra-l)S 
v + n — 1 

where v > controls the degree of moderation, and A > is the common 
k x k matrix toward which S is smoothed. In Section 4.1 we give the theo- 
retical reason for choosing this moderated variance-covariance matrix S and 
explain how we estimate v and A. Replacing S with S in the LiZ-statistic, 
our moderated X-R-statistic is 

(3.2) LR = 2(Z£ ax - Z£ ax ) = nlogf 1 + -^-d'ST 1 ^. 

When all the genes have an equal number of replicates n, equation (3.2) 
is a monotonic increasing function of ?7,d'S~ 1 d. We define the quadratic 
form nd'S _1 d = ||n 1//2 S _1 / 2 d|| 2 to be the moderated one-sample Hotelling 
T 2 -statistic. In the case of the null H : jjb = 0, S > 0, this is identical to the 
T 2 -statistic we derive in Section 4.1. The one-sample moderated Li2-statistic 
and the moderated Hotelling T 2 -statistic are hybrids of likelihood and Bayesian 
statistics, since S is estimated using the multivariate empirical Bayes pro- 
cedure we describe below. 

3.2. Unpaired two-sample problem. Similarly, in the unpaired two-sample 
case, the moderated Li?-statistic can be written as a function of the mod- 
erated two-sample Hotelling T 2 -statistic 

LR = (m + n) log ( 1 + ^l^d'S^d) , 

V m+n—2m+n J 

where d = Z — Y is the difference between sample averages and S = (m + 
n — 2) _1 ((m — l)Sz + (n — l)Sy) is the pooled sample variance-covariance 
matrix, and 

~_ (m + n -2)S + wV 

o — . 

m + n — 2 + v 

The term (m+n)~ 1 mnd'S~ 1 d is our moderated two-sample Hotelling T 2 -sta- 
tistic. We use the same approach to estimate S here as that for our two- 
sample multivariate empirical Bayes model described in Section 4.2. 
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4. The multivariate empirical Bayes model. 

4.1. One-sample or paired two-sample problem. 

4.1.1. Models and priors. The data Xi, . . . ,X n are multivariate normal 
with mean /x and covariance matrix S, denoted by Nk((i, S). We define an 
indicator random variable I to reflect the status of the gene, 

1, if if is true, 
0, if H is true. 



I 



We suppose that I has a Bernoulli distribution with success probability p, 

< p < 1 . Now we build up our multivariate hierarchical Bayesian model 
by first assigning independent and identical inverse Wishart priors to the 
gene-specific covariance matrices X: 

(4.1) S-Inv-Wishart^zM.)- 1 ), 

where v > and vK > are the degrees of freedom and scale matrix, respec- 
tively. Given S, we assign multivariate normal priors for the gene-specific 
mean /x for the two cases (I = 1) and (1 = 0): 

jU |S,/ = l~iV fe (0,7 ? - 1 S), 

fi\S, 1 = = 0, 

where r] > is a scale parameter. 

The posterior odds are the probability that the expected time course /x 
does not stay at (i.e., I = 1) over the probability that /x stays at (i.e., 

1 = 0), given the data Xi, . . . ,X n . Following [24] 's notation, we write 

P(J=l|Xi,...,X ra ) 

, P(/ = 0|X 1 ,...,X n ) 

4.2 

p P(Xi,...,X n / = l 



l-pP(X 1 ,...,X n |I = 0) 
The distribution of the data given I can be written as 

(4.3) P(X 1 ,...,X n |I)= fp(X|£,I)P(S|£,I)P(£|iV£. 

4.1.2. Multivariate joint distributions. Once the priors and the models 
are set, the joint distributions of the data can be determined given /. We 
omit the standard calculations leading to 

P(X 1 ,...,X„|/=1) 

r*((n + i/)/2) 



r fc ((n-l)/2)r fc (i//2) 



(4.4) 



Y. C. TAI AND T. P. SPEED 

X (n- l)(V2)fe(n-l) zy -(l/2)fcn 7r -(l/2)fc^-l + -lj-(l/2)fc 
| A |-(l/2)n| S |(l/2)(n-fe-2) 



life + {(n- 1 + Ti-^uA^XX! + (i/A/(n - l^-iSp/a)^) 

Thus, given J = 1, the probability density function of the data is a function 
of X and S only, which follows a Student-Siegel distribution [1]. Following 
[l]'s notation, this distribution is denoted by StSik{v 7 0, (n _1 +r/ _1 )A,n — 
l,(n - l)" 1 ^). For the case of / = 0, we get the same distribution with 
different parameters, namely, StSik(v, 0, n _1 A, n — l,(n-l) -1 i/A). 

4.1.3. MB-statistic and T 2 -statistic. Define our moderated gene-specific 
sample variance-covariance matrix S to be the inverse of the posterior mean 
of X -1 given S, 

(4.5) S = [E(^\S)]^= {n - 1)S + uA . 

y ' n-l+u 

The posterior odds O we defined earlier can be derived using the distribu- 
tions of the data given / and is 



O 



P 



v \(l/2)* 



1 — p \n + 7? . 

(4.6) 

n _l + u + T 2 X(l/2)(n+^) 



.n- l + v+{r]/{n + ifj)T 2 , 
where T 2 = t't and t is the moderated multivariate ^-statistic defined by 

(4.7) t = n 1 / 2 S" 1 / 2 X. 

Following the tradition in genetics, the log base 10 of O is called the LOD 
score. To distinguish it from the LOD score (also called the £>-statistic) in 
the univariate model of [20] and [24], the multivariate LOD score in this 
paper is called the MB-statistic, 

(4.8) MB = log 10 O. 

When all genes have the same number of replicates n, equation (4.8) is 
a monotonic increasing function of T 2 . This shows that the MB-statistic is 
equivalent to the T 2 -statistic when n is the same across genes, and therefore, 
one is encouraged to use the T 2 -statistic in this case since it does not require 
the estimation of r/ and leads to the same rankings as equation (4.8). We 
now derive the distribution for T 2 . 

By Gupta and Nagar [16], the Jacobian transformation from X to t is 
J(X^t) = |n- 1 /2s 1 /2|. 

Since equation (4.4) is a function of X and S only, 
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it is the joint probability density function for these two random variables. 
Substituting for X in terms of t in equation (4.4), and multiplying the 
resulting expression by J(X — > t), the joint probability density function for 
t and S given I = 1 is 



P(t,S|/ = l) 

= ,-(V^ r((n + y )/2) /n + gN _ 
r((n + i/-fe)/2) V n ) y 

(4.9) x 1 + ' (-5-W 

V n — 1 + z/ \ra + i] J 

1 



/2)fc 



/3 fc ((n-l)/2,i//2) 

| S |(l/2)(n-fe-2) 



|l/A/(Tl - l)|(V2)(n-l) | Ifc + (yX/(n _ l))-l S |(l/2)(n+v-l) " 

The above expression is factorized into parts involving S only and t only, 
proving that t and S are independent. It is apparent that t has a multivariate 
t distribution with degrees of freedom n + v — k, scale parameter n + v — 1, 
mean vector and covariance matrix r]~ l {n + irfjl^. This distribution is de- 
noted by t| I = 1 ~ t fe (n + i/ — k,n + v — 1, 0, ?y~ 1 (n + f?)Ifc) [16]. It is straight- 
forward to see that t\I = ~ t&(n + u — k,n + v — 1,0, I&). Given 1=1, 
S is distributed as a generalized type-II beta distribution with parameters 
(n — l)/2, v/2, scale matrix z/A/(ra — 1) and location matrix 0. The distri- 
bution is denoted by GBj} \(n — l)/2,u/2, uA/(n — 1), 0) [16]. The marginal 
distribution of S does not depend on I so that -P(S|/ = 0) = -P(S|I = 1). 
This distributional result is used to estimate the hyperparameter A. The 
distribution for T 2 under the null follows immediately. Under H, k~ lr t 2 
has an F distribution with degrees of freedom (fc,n + v — k); equivalently, 
(n + v — k)~ l {n + v— 1)T 2 has the Hotelling redistribution T 2 (k,n + u — 1). 

The T 2 -statistic is identical to the one-sample moderated Hotelling ^-sta- 
tistic in Section 3.1 with the same null hypothesis. 

For the easy extension to the above model, H : /i. = fi , XI > and K : fi ^ 
H , £ > 0, where fi Q is known, all the results above hold with X replaced by 
X-/x . 

4.1.4. Special cases. 

1. X = <7 2 Ifc. By constraining X! = <7 2 Ifc, we ignore the correlations among 
gene expression values at different times, and assume the variances at dif- 
ferent times are equal. Suppose that the prior for a 2 is 



a 2 ~ inv-gamma(^z/, ^u\ 2 ). 
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Define 

n 

^(n-ir'EP^-Xi) 2 ' 

1=1 

s| = (n - 1 + v)- l ((n - l)sf + u\ 2 ) 

and 

i,- =n 1/2 XjSj 1 , j = l,...,k. 

In this case, the posterior odds are equivalent to a product of k independent 
univariate odds, 

p ( _n_ \ {1/2)k A ( n-i + u + tj 2 \ (V2)(«+") 

1 - p \ra + r]J ^ Vn - 1 + 1/ + (^/(n + ?7))tj / 

(4.10) 

and the MB-statistic is equivalent to the sum of k univariate ^-statistics. 

2. n= 1. When n = 1, that is, when there is no replication at all, each 
gene has its own unknown variability. The moderated multivariate t-statistic 
becomes t = A~ 1 ^ 2 X. The posterior odds are obtained by plugging in n = 1 
in equation (4.6), and are found to be a function of X only. Since there is 
no replication, our hyperparameters must be assigned values, for example, 
from previous experiments. 

3. k = 1. When k = 1, that is, when there is only one time point, the 
alternative hypothesis states that there is differential expression at this single 
time point. Our multivariate model should and does reduce to the univariate 
model in [20] and [24]. 

4.1.5. Limiting cases. 

1. v — > oo. In this case, the gene-specific variance-covariance matrices 
are totally ignored. The moderated multivariate i-statistic above becomes 
too = n 1 / 2 A _1 / 2 X, and = t'^^. The posterior odds become 

I -p\n + r]J \2\n + r]J J 

2. v — > 0. In this case, there is no moderation at all. The posterior odds 
are just equation (4.6) with v replaced by 0. If n < k, then S _1//2 should be 
calculated by a ^-inverse. 
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3. v — > oo and 5] = a 2 Ik- Define iooj = n 1//2 A 1 X J -, j = l,...,k. The 
posterior odds become 




4. v — > and XI = o~ 2 Ik- In this case, the posterior odds are just equa- 
tion (4.10) with v replaced with 0. 

4.1.6. Hyperparameter estimation. We have shown that the MB-statistic 
for assessing whether or not a time course has mean depends on (k 2 + 
k + 6)/2 hyperparameters: u, A, rj and p. In practice, we need to estimate 
these hyperparameters, and plug in our estimates into the formulae for S, 
t, O, . . . and so on. Slightly abusing our notation, we will use the same 
symbols for these estimates, relying on context to make it clear whether we 
are assuming the hyperparameters to be known or not. In our multivariate 
model, many more hyperparameters need to be estimated, compared to the 
univariate models in [20] and [24], both of which have four hyperparameters. 
Closed form estimators for the hyperparameters in the univariate linear 
model setting are derived in [24], using the marginal sampling distributions 
of the statistic t and the sample variance s 2 , and are shown to be better 
than the simple estimators in [20] . Following [24] , the aim of this section is 
to derive estimators for the hyperparameters in our multivariate model. In 
general, the hyperparameter r\ associated with the case 7 = 1 is estimated 
based on only a small subset of genes, while v and A are estimated using 
the whole gene set. Instead of estimating the proportion of differentially 
expressed genes p, we plug in a user-defined value, since the choice of p does 
not affect the rankings of genes based on the MB-statistic. 

EB estimation of u and A. The hyperparameter v determines the degree 
of smoothing between S and A. The method we use to estimate v builds on 
that used to estimate do in Section 6.2 in [24]. However, unlike do in [24], v 
is associated with the k x k matrix XI. Therefore, a method appropriate to 
this multivariate framework is needed. Let Uj be the estimated prior degrees 
of freedom based on the jth diagonal elements of the gene-specific sample 
variance-covariance matrices (i.e., the replicate variances for the jth time 
point from the whole gene set) using the method proposed in Section 6.2 in 
[24]. Our estimation of v is based on the following two-step strategy. As the 
first step, set v as v = max(mean(i>j), k + 6),j = 1, . . . ,k. This estimated u 
is used to estimate A. Once A is estimated, v is reset to be v = mean(i)j). 
In practice, one can even just plug in a user-defined value vq which gives 
the desired amount of smoothing. In such the first step sets u = 
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max(fo, k + 6). This v is used to estimate A. After A is estimated, v can be 
reset to the user-defined value vq. 

Our estimate of A comes after the first step in the estimation of v . We 
showed that, under our model, S follows the generalized type-II beta distri- 
bution with expectation [y — k — l)~ 1 z^A. By the weak law of large numbers, 
S converges in probability to {y — k— l) -1 z^A. We can thus estimate A by 
z)~ 1 (z) — k — 1)S. If v — ► oo, then A is estimated by S. The above estimates 
give quite satisfactory results on real data. A theoretical analysis of the es- 
timation of our hyperparameters will be given later. For the moment, we 
content ourselves with obtaining reasonable estimates. 

EB estimation ofrj. The hyperparameter r] is related to the moderated 
multivariate t-statistic t of nonzero genes. The method we use to estimate 
77 builds on that of estimating vq in [24], except that we now need to deal 
with the multivariate case. Let tj be the jth element of the moderated 
multivariate i-statistic t, j = 1, . . . , k. As in Section 6.3 in [24], each tj gives 
an estimate of rj, call it rjj, based on the top p/2 portion of genes with the 
largest \tj\. We set 17 to be the mean of the r)j. 

4.2. Unpaired two-sample problem. Suppose there are two independent 
biological conditions Z and Y with sample sizes m and n, respectively. 
We can also derive the MB-statistic for testing the null H : = Yl z = 
Sy > 0. The null hypothesis turns out to be the same as that in Section 4.1: 
H : ju = 0, 5] > 0, if we write fi = fi z — n Y and S = Yl z = Sy. That is, we 
solve this two-sample problem using the one-sample approach in Section 4. 1 . 
We denote the m i.i.d. time course vectors for biological condition Z by 
Zi, . . . , Z m , each from a multivariate normal distribution with mean ^ z and 
variance-covariance matrix S. Similarly, those for biological condition Y are 
denoted by Yi, . . . , Y n , each with mean fi Y an d variance-covariance matrix 
5]. Since the null hypothesis here is identical to that in Section 4.1, the 
priors for \x and £ are exactly the same as those in Section 4.1, and we omit 
the details here. In a later paper we will attack this problem by assigning 
independent priors for fj, Y and fi z separately. 

All the results follow immediately. The moderated multivariate f-statistic 
t here is defined as equation (4.7) with n replaced by (m _1 + n~ l )~ l and X 
replaced by Z — Y. S here is the same as that defined in Section 3.2. The 
posterior odds O against the null hypothesis that the expected time courses 
are the same are 




( 



m + n-2 + v + T 2 



(l/2)(m+n+u-l) 



X 



m + n — 2 + v + {{m 1 +n 



l )/(m- 1 +n- 1 + n- l ))T 2 
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The log base 10 of O is our two-sample MB-statistic. Again, when all genes 
have the same sample sizes m and n, the two-sample MB-statistic is equiv- 
alent to the T 2 = t't. Under H, k~ l T 2 has an F distribution with degrees 
of freedom (k, m + n + v — k — 1); equivalently, [m + n + v — k — 1) _1 (m + 
n + v — 2)T 2 has the Hotelling T 2 distribution T 2 (k, m + n + v — 2). 

The MB-statistic described in this section involves hyperparameters v, 
A, rj and p. The estimation procedures for these hyperparameters are very 
similar to those in Section 4.1, except that we have to use the gene-specific 
pooled sample variance-covariance matrices when estimating v and A, and 
use the x 1 moderated multivariate i-statistic t defined here to estimate 
r]. We omit the details here. 

It should be noted that the MB-statistic derived in this section has a 
slightly different definition: instead of using all the data we observe, we 
only use the difference in sample averages and the pooled sample variance- 
covariance matrix. The T 2 -statistic here is identical to the moderated two- 
sample Hotelling T 2 -statistic in Section 3.2. 

4.3. One-sample problem of constancy. In this section we derive the pos- 
terior odds against the null that a gene's mean expression level stays constant 
over time. We obtain a closed form solution similar to that in the preceding 
sections, but only under a constraint on the variance-covariance matrix E. 

4.3.1. Transformation. For each gene, let / be the indicator variable 
defined in Section 4.1. Let P = be the k x k projection matrix onto 

the rank 1 space of constant vectors, where ll = (1, . . . , 1) is a k x 1 vector 
of Is. Let P c = Ifc — P be the projection onto the orthogonal complement of 
-R(P). We can write any vector \x 6 R k as /x = P/x + P c /x, and in the case 
1 = 0, the second term P c /x vanishes. As in Section 4.1, we build up our 
multivariate model by first assigning independent inverse Wishart priors to 
the gene-specific covariance matrices 'S; see equation (4.1). Given 5], we next 
assign multivariate normal priors to the gene-specific mean parameters /x for 
the case of nonconstant (/ = 1) and constant genes (1 = 0), respectively: 

, , f //|I],/=l-Af(0,r- 1 PSP + K - 1 P c SP c ), 

1 ' ' \ /i|E,/ = 0~JV(0, t _1 PEP). 

Given X and 1 = 0, the covariance matrix PEP guarantees that /z is a 
constant vector, while when 1=1, the extra component P C SP C adds fur- 
ther variance to n so that it becomes a nonconstant vector. Again, in or- 
der to obtain the full expression for the posterior odds O, we need to de- 
rive P(Xi, . . . ,X n |7) using equation (4.3). To get a closed-form expression 
for the posterior odds, we find it necessary to make an additional assump- 
tion, namely, that PE = EP. With this assumption, given E and / = 0, X 
is a multivariate normal distribution with mean and covariance matrix 



14 



Y. C. TAI AND T. P. SPEED 



(n _1 E + t _1 SP). Similarly, given E and 1=1, X is a multivariate nor- 
mal distribution with mean and the covariance matrix (n -1 E + r _1 SP + 

K^EP ). 

For the rest of this section, unless stated otherwise, we assume PS = SP, 
and we make use of the following lemma, whose proof is omitted. 

Lemma 4.1. Suppose T is any k x k nonsingular matrix whose first row 
is constant c and the remaining rows have row sums equal to 0. Write T = 
(Tq,T^)', where Tq is the first row of T, and Ti is the remainder. Then, 
for any E > satisfying PS = SP, TST' = S is a k x k block diagonal 
matrix with the scalar a 2 > as the first block and (k — 1) X (k — 1) matrix 
Si > as i/ie second block: that is, 

As the first example, let T be the Helmert matrix, where the jith. element 
of T is defined as 

' tji = 1/y/k, for j = 1, i = 1, . . . , k, 

i <j< = Wj(j-1 ), for 2<j<fc,l<i<J-l, 

<;< = -0'- Wj'O'-l), ^ 2<j<k,i = j, 

>tji = 0, for 2<j <k-l,j + l<i<k. 

T can also be the following matrix, where the jith element of T is defined 

as 

1, for j = l,i = l,...,k, 
1, for 2 < j < k, i = j — 1, 
-1, for 2< j < k,i=j, 
0, otherwise. 

For our multivariate empirical Bayes model in this section, we use the 
Helmert matrix T to proceed with our calculations. The results can be 
applied to other T immediately. 

4.3.2. Models and priors. Here T is partitioned into its first row To 
(1 x k) and its last k — 1 rows Ti {{k — 1) x k). Since Xi, . . . ,X n are i.i.d. 
Nk(fi,Jj), the transformed random vectors TXj are also multivariate nor- 
mally distributed with mean T/z and covariance matrix S. By Lemma 4.1, 
the matrix S is a block diagonal matrix with a 2 as the first block, and Si 
as the second block. Defining X{ = k" 1 Ylj=i Xiji then \/kxi and the random 
vectors TiXj are independent and normally distributed, with distributions 

Vkxi\T fi,a 2 ~iV(T /x,a 2 ), 

TiX^T^i, Si Ei). 




REPLICATED MICROARRAY TIME COURSE DATA 



15 



This transformation allows us to separate the gene expression changes into 
constant and nonconstant changes. 

As we have seen in Section 4.1, the joint distributions of data given / can 
be fully described using the sufficient statistics x, TiX, s 2 and Si, where 
x = n- 1 J2 7 Li^ TiX = n" 1 Er=iTiX 4 , s 2 = (n - l^E^a* -J) 2 and 
Si = (n - I)" 1 E?=i(TiXi - TiX)(TiXj - TiX)'. The prior for S is first 
set through the independent priors for a 2 and Si. We suppose that a 2 and 
Si are independently distributed, with an inverse gamma distribution with 
shape parameter £/2 and scale parameter £A 2 /2, and an inverse Wishart 
distribution with degrees of freedom v and scale matrix z^Ai, respectively, 
that is, 



(4.12) 



(4.13) 



a 2 ~ inv-gamma(|£, ^£A 2 ), 
Si ~ inv-Wishart^i/Ai)- 1 ). 

The prior for Tfi has four parts. We assign independent priors to To/x and 
Ti/^t separately for the cases 1 = 1 and 1 = 0. For the case 1=1, priors are 

T n\a 2 ,I= l~N(6,n- 1 dl), 
Ti/x|Si,/=l~iV(0,7 ? - 1 Si), 

where 9 > is the mean, and k > and rj > are scale parameters. When 
1 = 0, Ti/x = with probability 1. Thus, the priors in this case are 

T /i|cr 2 ,/ = 0~ N{6,n~ l d 2 ), 
Ti/i|Si, 1 = = 0. 

It is reasonable to assume P(To/x|<r 2 , / = 0) = P(Tofi\a 2 , 1 = 1) for large 
genome- wide arrays since there is no obvious reason why the expected grand 
mean of the expression levels for nonconstant genes should differ from that 
of constant genes. For two-color comparative microarray experiments, it is 
also reasonable to assume 6 = 0. 



(4.14) 



4.3.3. Multivariate joint distributions. The joint distributions can be de- 
rived quite readily using a previously established formula, and so we omit 
the details. Given 1=1, TiX and Si follow the Student-Siegel distribu- 
tion StSik-i{v,Q, (n _1 + r/ _1 )Ai,n — 1, (n — l) _1 z^Ai). Similarly, the joint 
distribution of TiX and Si given / = is StSik-\(v, 0, n _1 Ai, n — 1, (n — 
l)-^Ai). 

4.3.4. MB -statistic and T 2 -statistic. The posterior odds against the null 
that a gene's mean expression level stays constant over time are equation 
(4.6) in Section 4.1 with k replaced by k — 1, t expressed by equation (4.7) 
with S replaced by Si and X replaced by TiX. Si is just equation (4.5) with 
S replaced by Si and A replaced by Ai. As in Section 4.1, the MB-statistic 
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is a monotonic increasing function of T 2 = t't when all genes have the same 
sample size n. 

Under H, (k — 1) _1 T 2 has an F distribution with degrees of freedom 
{k — l,n + u — k + 1), or, equivalently, (n + v — k + l) _1 (n + v — 1)T 2 has 
a Hotelling T 2 -distribution T 2 (k — l,n + v — 1). The hyperparameter esti- 
mation procedures here are similar to those described in Section 4.1, except 
that all the estimations are performed based on transformed data. 

5. Case study. In this section we illustrate our results with a paired two- 
sample problem, using the Arabidopsis thaliana dataset in [35]. Here we only 
give a very brief description of the data and the results. We refer the reader 
to [35] and to Chapter 5 of [28] for more thorough discussions. 

A. thaliana wildtype (Columbia) and icsl-2 null mutant plants were evenly 
positioned, intermixed and grown in growth chambers under controlled con- 
ditions. When the plants were four weeks old, they were infected with a 
moderately heavy innoculum of the powdery mildew G. orontii. Each pair 
of mRNA samples from wildtype and mutant plants was harvested and col- 
lected at six time points post-infection. Plants could not be resampled, so 
mRNA samples at one time point were from different plants than those 
of any other time point. We report here on the analysis of three replicate 
experiments under similar environmental conditions which contribute four 
biological replicates: one from the first and third experiments (1-3, 3-1) 
and two from the second experiment (2—1, 2-2). These mRNA samples were 
hybridized onto Affymetrix Arabidopsis ATH1 GeneChips, yielding 22,810 
probesets and 48 arrays in our analysis. The array preprocessing were done 
using the Robust Multi-array Analysis (RMA) algorithm described in [5, 18] 
which is implemented in the Bioconductor package af fy. 

This study is longitudinal if we treat experiments as units, while it is cross- 
sectional if we treat plants as units. We believe it is worthwhile to treat it as 
a paired longitudinal study, since samples within the same experiment are 
more similar than those from different experiments. We thus have a paired 
two-sample problem, with the genes of interest being those whose wildtype 
and mutant temporal profiles are different. We subtracted the log 2 intensities 
of the wildtype from those of the paired mutant at each time within each 
replicate, yielding the log 2 ratios of mutant relative to wildtype. Since the 
number of time course replicates is the same [n = 4) across genes for this 
dataset, we used the T 2 -statistic instead of the MB-statistic to rank genes, 
so that we did not have to estimate the hyperparameter rj. 

For comparison purposes, we fitted a linear model to the log ratios for 
each gene, with time and replicate effects, and calculated the F-statistic for 
the time effect. 



REPLICATED MICROARRAY TIME COURSE DATA 



17 



Table 1 

Spearman rank correlations between T 2 with different v 
and the estimated v for all and the top 859 genes. The 
percent moderation is defined by (y + n — 1)~ v x 100 



% moderation 


Correlation (all) 


Correlation (top 859) 


97 = 100) 


0.97 


0.90 


80 = 12) 


0.99 


0.98 


40 [v = 2) 


0.99 


0.98 


25 \v = 1) 


0.98 


0.96 


o = o.oi) 


0.93 


0.90 



5.1. Results. The extent of moderation from i> = 5 was 63%. The left 
panel of Figure 1 presents three genes falling into different ranges of ranks 
(rank = 1, 175, 859) by T 2 , while the ones on the right panel have the same 
ranks by F. The gene ranked most highly by T 2 exhibits much greater 
differences between the wildtype and mutant temporal profiles than the one 
ranked most highly by F. The magnitude of the difference, as measured 
by T 2 , decreases as the rank goes down. The gene ranked 1 by T 2 is well 
known: pathogenesis-related protein 1 (PR1). Other known pathogenesis- 
related genes also ranked highly by T 2 and were less highly ranked by the 
F-statistic, as detailed in [35]. 

To investigate the effect of the amount of moderation on gene ranking, 
we kept A fixed, and re-calculated the T 2 -statistic with several different 
v values. We then computed the Spearman rank correlation between the 
different sets of T 2 -statistics for all genes and for the top 859 genes only. 
Table 1 gives the results. The correlations are lower in the two extremes. All 
of these sets have the same number one gene. This comparison shows that 
the gene ranks are reasonably stable when the extent of moderation varies 
within a certain window, and that moderation seems to have more effect on 
the top genes relative to the whole gene set. 

6. Simulation study. 

6.1. Method. In this section we report on a small simulation study for 
the null hypothesis H : /x = hq\, S > based on an actual example we have 
met. We simulate 100 data sets, each with 20,000 genes. The genes are 
simulated independently, which we regard as an assumption that makes sense 
to compare methods, but it should be kept in mind that gene expression 
measures in real data can be quite dependent. In each simulated data set, 
400 out of the 20,000 genes are assigned to be nonconstant. That is, p = 0.02. 
Each gene is simulated with three independent replicates (n = 3) and eight 
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Fig. 1. Genes of ranks 1, 175, 859 by the paired two-sample T 2 -statistic in Section 4.1 
(left panel) and the F-statistic (right panel) . The temporal difference between wildtype and 
mutant decreases as the rank by T 2 goes down. The genes on the left panel all show larger 
differences between the wildtype and mutant than the corresponding ones on the right panel. 
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time points {k 



The other hyperparameters are the following: v = 13, 



£ = 3, A = 0.3, 9 = (two-color experiments), k = 0.02, rj = 0.08, and 



( 14.69 
0.57 
0.99 
0.40 
0.55 
0.51 

V-0.23 



0.57 
15.36 
1.22 
0.84 
1.19 
0.91 
0.86 



0.99 
1.22 
14.41 
2.47 
1.81 
1.51 
1.07 



0.40 
0.84 
2.47 
17.05 
2.40 
2.32 
1.33 



0.55 

1.19 

1.81 

2.40 

15.63 

3.31 

2.75 



0.51 

0.91 

1.51 

2.32 

3.31 

13.38 

3.15 



-0.23\ 

0.86 

1.07 

1.33 

2.75 

3.15 
12.90 / 



x 10" 



The correlation matrix of A is 



0.04 


0.07 


0.03 


0.04 


0.04 


-0.02 




1 


0.08 


0.05 


0.08 


0.06 


0.06 


0.08 


1 


0.16 


0.12 


0.11 


0.08 




0.05 


0.16 


1 


0.15 


0.15 


0.09 




0.08 


0.12 


0.15 


1 


0.23 


0.20 




0.06 


0.11 


0.15 


0.23 


1 


0.24 




0.06 


0.08 


0.09 


0.19 


0.24 


1 


/ 



1 

0.04 
0.07 
0.03 
0.04 
0.04 
V-0.02 

and we see clear evidence of serial correlation. Note that, in the real world, 
although often the case, the correlation does not always decrease with time 
lag. The statistics compared in our study are the following: 

(1) MB-statistic, or equivalently, the T 2 -statistic; 

(2) MB-statistic using first differences: take the differences in gene ex- 
pression values at consecutive time points within replicates, and use them 
to test the null hypothesis H : /i. = 0, X > (Section 4.1), where fi is the 
mean of the differences; 

(3) MB-statistic in the special case 5] = <J 2 Ik'> 

(4) MB-statistic in the limiting case v — ► oo; 

(5) MB-statistic in the limiting case v — > 0; 

(6) ordinary F-statistic from an ANOVA model with time and replicate 
effects; 

(7) partly-moderated F-statistic proposed in [24] from an ANOVA model 
with time and replicate effects; 

(8) one-sample moderated Hotelling T 2 -statistic ||n 1 / 2 S~ 1 / 2 d|| 2 derived 
in Section 3, equivalently, the moderated BB-statistic, where the degree of 
moderation and the common matrix toward which each sample covariance 
matrix moves is estimated by the method given in Section 4.1; 

(9) the variance across time course replicates (nk - l)" 1 £?=l EjUPQj - 



Here each of the nine statistics incorporates either none (e.g., variance) or 
one (ordinary F-statistic) or more of the following: moderation, correlation 
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Table 2 

The means and standard deviations (SD) of the diagonal 
elements of the estimated Ai 



Hyperparameters 


True value X 10 3 


Mean X 10 3 


SD xlO 3 


A? 


14.69 


14.71 


0.16 




15.36 


15.37 


0.17 


\2 

A'i 


14.41 


14.43 


0.15 


\2 
A 4 


17.05 


17.04 


0.19 


A B 


15.63 


15.63 


0.15 


A§ 


13.38 


13.40 


0.15 


A? 


12.90 


12.92 


0.17 



structure and replicate variances, and thus can be used to show the impor- 
tance of the above properties. It is not appropriate to set the prior degrees 
of freedom v to be a very small number, since we have the constraint that 
v > k — 1. We choose v to be k + 5 = 13 because it simulates more stable 
S's across genes. 

6.2. Results. Table 2 compares the means and standard deviations of 
the hyperparameter estimates of the diagonal elements of Ai (A|), j = 
1, . . . , k — 1, with their true values. The mean estimate of Ai is very close to 
the true Ai, and the standard deviations are very small. The hyperparameter 
j) is always under-estimated (mean = 0.026, SD = 0.002), which agrees with 
Section 8 in [24], where vq was usually over-estimated. The hyperparame- 
ter v is also always under-estimated (mean = 7.0, SD = 0.2). In Section 5 
we observed that the amount of moderation v does not greatly affect gene 
ranking except at the two extremes. One can even choose a user-defined v 
which gives reasonable results. Although not well estimated, r\ only affects 
the rankings when the number of replicates is different across genes. How- 
ever, this does not happen often in the real world. Even when that happens, 
the effect is very small. To investigate the effects of rj on gene rankings, we 
tried a couple of r/ values from different ranges, while keeping the remaining 
hyperparameters fixed, and calculated the MB-statistics. The rank correla- 
tions between rankings of the MB-statistics with the user-defined ry's and 
the estimated rj for one set of simulated data are the following: 0.91, 0.94, 
0.99, 0.99, 0.99 for rj=2, 1, 0.08 (true value), 0.05, 0.001, respectively. 

To examine the relationship between the T 2 -statistic and the true de- 
viation from constancy, the log 10 transformed T 2 -statistic from one simu- 
lated dataset is plotted against the Mahalanobis distance between the ex- 
pected time course vector n and its projection onto the rank 1 constant 
space p. = P/^ (Figure 2). The squared Mahalanobis distance is defined by 
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d(/x, p,) 2 = (fi — — pi). Figure 2 clearly shows that log 10 T 2 is pos- 

itively correlated with d(fj,,fj,), and most of the 400 true nonconstant genes 
achieve higher T 2 -statistics than the constant genes. 

Figure 3 plots the average number of false positives against average num- 
ber of false negatives at different cutoffs. By different cutoffs, we mean choos- 
ing the top x genes and calculating the numbers of false positives and false 
negatives, where x varies across the integers from 400 to 800. The lines 
in Figure 3 from left to right represent the following: MB-statistic (T 2 ), 
MB-statistic with first differences (indistinguishable from the MB-statistic), 
one-sample moderated Hotelling T 2 -statistic (indistinguishable from the 
MB-statistic), MB-statistic with 5] = <r 2 Ifc, MB-statistic with v — > oo, partly- 
moderated F-statistic [24], ordinary F-statistic, MB-statistic with u — > 
and variance. The MB-statistic (T 2 ) attains almost the same numbers of 
false positives and false negatives as MB with first differences and the one- 
sample moderated Hotelling T 2 -statistic. The effectiveness of moderation is 
demonstrated by comparing the lines for the MB-statistic, the MB-statistic 
in the limiting case v — > oo and the MB-statistic in the limiting case that 
v — ► 0. Both of these limiting cases lead to higher aggregate false positives 
and false negatives. This result supports the view stated in [29] that moder- 
ation is useful. In particular, the case v — > (no moderation at all) produces 




20 40 60 80 



Fig. 2. The log 10 T 2 statistic versus the true deviation from constancy d(fJL,jl) for one 
simulated dataset. Here 1 denotes nonconstant, and o constant genes. 
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Fig. 3. Average number of false positives versus number of false negatives of all the nine 
statistics. The subplot presents the curves for the best seven statistics. 



much higher numbers of false positives and false negatives. This is likely 
due to the poor estimation of sample variance-covariance matrices with a 
small number of replicates. Indeed, the ordinary unmoderated F-statistic 
which ignores the correlation structure achieves smaller numbers of false 
positives and false negatives than the unmoderated MB-statistic. A similar 
situation also arises in the microarray discrimination context; see Section 7 
of [12]. The partly-moderated F-statistic [24] which ignores the dependence 
among times behaves like the MB-statistic in the special case S = cr 2 Ifc. 
Moreover, it achieves fewer false positives and false negatives than the or- 
dinary F-statistic. Figure 3 also demonstrates the importance of incorpo- 
rating the correlation structure among time points. The MB-statistic (T 2 ) 
and the one-sample moderated Hotelling T 2 -statistic perform better than 
the partly-moderated F-statistic in [24] and the ordinary F-statistic; the 
former incorporates the correlation structure among time points, whereas 
the latter does not. However, we observe that the amount of moderation 
given by the partly-moderated F-statistic in [24] is usually much less than 
that given by the MB-statistic. When there are a large number of residual 
degrees of freedom from the linear model, the partly-moderated F-statistic 
[24] behaves very much like the ordinary F-statistic. This suggests that 
the lower number of false positives and number of false negatives from the 
MB-statistic than the partly-moderated F-statistic in [24] involve both the 
incorporation of correlation structures and the amounts of moderation. As 
expected, the simple variance statistic across replicates, which totally ig- 
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nores the replicate variances, performs the worst. This demonstrates the 
importance of incorporating the replicate variances into any statistic. 

7. Discussion. In this paper we introduced the MB- and T 2 -statistics 
for one- and two-sample longitudinal replicated developmental microarray 
time course experiments. Our main focus was the one-sample or paired two- 
sample problem with the null hypothesis H : fi = 0, S > 0. This MB-statistic 
can be used when there are two biological conditions and the samples are 
paired across conditions, and it is shown to perform better than the classi- 
cal .F-statistic on a problem briefly described in Section 5. In addition, we 
also derive the MB-statistics and T 2 -statistics for the two-sample problem 
with the null H : \x z = fj, Y , = = 5] > 0, and the one-sample problem 
with the null H : fi = (jlqI, £ > using similar approaches. The latter situa- 
tion requires a slight assumption on £ in order to get a simple closed-form 
solution for the posterior odds against the null. All the MB-statistics and 
T 2 -statistics incorporate the correlation structure, replication and modera- 
tion. The moderated versions of some standard likelihood-ratio test statis- 
tics are also described. When all genes have the same sample size(s), our 
T 2 -statistics are not only equivalent to the MB-statistics, but also are iden- 
tical to their corresponding moderated Hotelling T 2 -statistics, apart from 
the one-sample problem with the null H : /x = /jlqI, £ > 0, where there is an 
additional constraint on X. In this case the T 2 -statistic performed as well 
as the moderated Hotelling T 2 -statistic in our simulation study, and also on 
several real datasets we have encountered. We have shown in the simulation 
study that, with this null, the MB-statistic (T 2 ), the MB-statistic using first 
differences and the one-sample moderated Hotelling T 2 -statistic perform 
best among all the nine statistics compared. This is not entirely surprising 
given that we simulated data under our model, but the comparisons are still 
informative. In practice, we regard the MB-statistic, the MB-statistic with 
first differences and the moderated Xi2-statistic as performing equally well, 
giving very similar (or identical) results. However, to use the /^-statistic 
(the moderated Hotelling T 2 -statistic), we still need to insert moderated 
sample variance-covariance matrices, and these come from our multivariate 
empirical Bayes framework. In other words, our models provide a natural 
way to estimate the gene-specific moderated variance-covariance matrices 
(Sections 4.1-4.3), while the likelihood-based approach alone does not. 

The assumption of PS = SP with the null H: fj, = /iol> S > allows the 
mathematical calculations in Section 4.3, and leads to a closed-form for- 
mula for the MB-statistic. One question which naturally arises to be the 
impact of this constraint on the rankings of genes. From the practical point 
of view, the impact of this constraint on gene rankings is very slight. The 
rank correlations between the one-sample MB-statistic with the commuting 
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assumption and the moderated B-R-statistic from likelihood-based approach 
without the constraint, from the actual examples we have met, are typically 
very high (over 0.99). The rank correlations from our simulated data are 
also over 0.99. On the other hand, using the MB-statistic with first differ- 
ences produces very similar or even identical results to the MB-statistic in 
Section 4.3. Indeed, instead of using the Helmert matrix, if we choose T to 
be the second transformation matrix of Section 4.3, we get identical results. 
Even so, we do not consider the first differences approach to be the solution 
to this commuting constraint, since the inference drawn is based on reduced, 
not the original data. In other words, the null hypotheses are not equivalent. 
The likelihood-based approach with moderation and the first differences ap- 
proach support the fact that this constraint does not have much effect on 
the results. In addition, the former is a good way to avoid this assumption, 
since it performs as well as our MB-statistic in Section 4.3. 

The statistics proposed in this paper are for one- and two-sample longi- 
tudinal data. We should be aware that many experiments in the real world 
exhibit some features from both longitudinal and cross-sectional experiments 
(e.g., Section 5). 

One thing we plan to investigate in the future is the effect of assuming 
the same variance-covariance matrix S for both 1 = 1 and 1 = 0. Another 
issue which interests us is the effect of assuming the same S across biological 
conditions in the unpaired two-sample model in Section 4.2. The proposed 
methods may be extended in several ways, for example, identifying genes of 
some specific pattern, rather than any pattern. The statistics for a longitudi- 
nal multi-sample problem when there are at least three biological conditions 
and genes of interest are those with different temporal profiles across con- 
ditions derived in [30]. The corresponding statistics for cross-sectional data 
are also presented in [31]. 

The proposed methods focus on gene ranking, but not assessing the signifi- 
cance using p-values. However, if desired, we believe that generating p- values 
from a bootstrap analysis should be successful in this context. 

We constructed our models using conjugate priors for the multivariate 
normal likelihoods, so that we got simple closed- form solutions for the pos- 
teriors odds. Finding a closed-form statistic when the priors on fj, and S are 
independent seems to be an open and probably hard problem; that problem 
probably needs to be dealt with using MCMC. 
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